
clear;
close all;
%% select channel
P="D:\APPdocu\python\UWA_GAN\data";
GAN_Path=P;
flag=1; % 1: NOF1; 2: NSC1
channel1='Generated_data\NOF1_256';
channel2='Generated_data\NCS1_256'; 
filepath1 = fullfile(GAN_Path, channel1, 'mat');
filepath2 = fullfile(GAN_Path, channel2, 'mat');
% 
% if flag==1
%    fs_tau=16000;
%    fs_t=fs_tau/510/4;
% else
%    fs_tau=16000;
%    fs_t=fs_tau/2048;
% end

% temp=22;
% filename1 = fullfile(filepath1, ['mat_' sprintf('%03d', temp) '.mat']);
% filename2 = fullfile(filepath2, ['mat_' sprintf('%03d', temp) '.mat']);
% 
% var_struct1=load(filename1);
% name_cell=fieldnames(var_struct1);
% mat1=getfield(var_struct1,char(name_cell));
% 
% 
% var_struct1=load(filename2);
% name_cell=fieldnames(var_struct1);
% mat2=getfield(var_struct1,char(name_cell));
% 
% 
% h1_real=mat1(:,:,1);
% h1_img=mat1(:,:,2);
% h1=h1_real+1i*h1_img;
% 
% h2_real=mat2(:,:,1);
% h2_img=mat2(:,:,2);
% h2=h2_real+1i*h2_img;


%% 定义信道参数
%--------------------------------------------------------------------------

M_mod = 4;                         % 4QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 
fc = 12500;
c = 1500;                                                     % speed of light (m/s)
N = ofdmSym;

Fn = dftmtx(N);       % Generate the DFT matrix
Fn = Fn./norm(Fn);    % normalize the DFT matrix

Fn2 = dftmtx(M);       % Generate the DFT matrix
Fn2 = Fn2./norm(Fn2);    % normalize the DFT matrix

Md = M-Mzp;

%% 时延降采样且时间升采样
% [~,G1]=channel_acc(h1,M,N,1);
% [~,G2]=channel_acc(h2,M,N,2);
%% 发送信号
data_grid = zeros(M,N);
data_grid(1:Md,:) = 1;
N_syms_perfram = sum(sum(data_grid));

delta_f = scs;
T = 1/delta_f;
% length_CP=cpLen;

% max_speed=20;  % km/hr
% [chan_coef,delay_taps,Doppler_taps,taps]=Gen_delay_Doppler_channel_parameters(N,M,fc,delta_f,T,max_speed,7);
% L_set=unique(delay_taps);
% 
% [G,gs]=Gen_time_domain_channel(N,M,taps,delay_taps,Doppler_taps,chan_coef);


% [~,K_ml]=Gen_DD_channel_kml(N,M,L_set,gs); 
% [H_tf]=Gen_time_frequency_channel(N,M,gs,L_set);
             
%%                       OTFS BER Calculation

     data_na= randi([0,1],N_syms_perfram*k,1);
     data=qammod(reshape(data_na,k,Md*N), M_mod,'gray','InputType','bit');
     X = Gen_2D_data_grid(N,M,data,data_grid);
    
    
    %% OTFS modulation%%%%
        X1_tilda=X*Fn';
        s1= reshape(X1_tilda,N*M,1);
    
     %% channel output%%%%%

% %   r1=zeros(N*M,1);
%     noise1= sqrt(0.2/2)*(randn(size(s1)) + 1i*randn(size(s1)));
%     l_max=max(delay_taps);
%     for q=1:N*M
%         for l=(L_set+1)
%             if(q>=l)
%                 r1(q)=r1(q)+gs(l,q)*s1(q-l+1);
%             end
%         end
%     end
 
 
     %% OTFS demodulation%%%%
%     Y1_tilda=reshape(r1,M,N);
%     Y1 = Y1_tilda*Fn;
%     y1=reshape(Y1.',N*M,1);
%% detect
   eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
    SNR_dB = -5:10:25;
    SNR = 10.^(SNR_dB/10);
    noise_var_sqrt = sqrt(1./SNR);
    sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;     

    errorRate1 = comm.ErrorRate('ResetInputPort',true);
    errorRate2 = comm.ErrorRate('ResetInputPort',true);

    berOTFS1 = zeros(length(SNR),3); 
    berOTFS2 = zeros(length(SNR),3);

   
    for temp = 8
        disp(temp)
        filename1 = fullfile(filepath1, ['mat_' sprintf('%03d', temp) '.mat']);
        filename2 = fullfile(filepath2, ['mat_' sprintf('%03d', temp) '.mat']);
        
        var_struct1=load(filename1);
        name_cell=fieldnames(var_struct1);
        mat1=getfield(var_struct1,char(name_cell));
        
        
        var_struct1=load(filename2);
        name_cell=fieldnames(var_struct1);
        mat2=getfield(var_struct1,char(name_cell));
        
        
        h1_real=mat1(:,:,1);
        h1_img=mat1(:,:,2);
        h1=h1_real+1i*h1_img;
        
        h2_real=mat2(:,:,1);
        h2_img=mat2(:,:,2);
        h2=h2_real+1i*h2_img;
        [~,G1]=channel_acc(h1,M,N,1);
        [~,G2]=channel_acc(h2,M,N,2);
        for m = 1:length(SNR)
            noise= sqrt(sigma_2(m)/2)*(randn(size(s1)) + 1i*randn(size(s1)));
    
            r1_n=G1*s1+noise;
            r2_n=G2*s1+noise;
    
         
    
    %         figure(3)
    %         plot(0:1/B:(M*N-1)/B,abs(r1),'y',0:1/B:(M*N-1)/B,abs(r2),'r');
    %         legend('NOC1','NCS1')
    % %         xlim([0,0.4]);
    %         title('输出信号')
    
            [est_LMMSE_r1_n,~] = LMMSE_detector_OTFS(N,M,M_mod,sigma_2(m),data_grid,r1_n,G1); 
            [est_LMMSE_r2_n,~] = LMMSE_detector_OTFS(N,M,M_mod,sigma_2(m),data_grid,r2_n,G2); 
            berOTFS1(m,:)=berOTFS1(m,:)+errorRate1(data_na,est_LMMSE_r1_n,1).';
            berOTFS2(m,:)=berOTFS2(m,:)+errorRate1(data_na,est_LMMSE_r2_n,1).';
            
        end
      
    end

%% figure
berOTFS1(~berOTFS1)=1e-12;
berOTFS2(~berOTFS2)=1e-12;


figure
berOTFS1(:,1)=berOTFS1(:,2)./berOTFS1(:,3);
berOTFS2(:,1)=berOTFS2(:,2)./berOTFS2(:,3);

berOTFS1(~berOTFS1)=1e-12;
berOTFS2(~berOTFS2)=1e-12;

semilogy(SNR_dB,berOTFS1(:,1),'-r');             %Plot simulated BER w/ OTFS
hold on;
semilogy(SNR_dB,berOTFS2(:,1),'--b');
% semilogy(SNR,berMChannel(:,1),'--b');             %Plot simulated BER w/ C-OTFS


ylabel('BER');
xlabel('SNR/dB');
legend( 'NOF1',"NCS1");
title("输出误码率")
grid on;
hold off;
